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ABSTRACT 

In this paper we review the development of the shock-capturing methodology, 
paying special attention to the increasing nonlinearity in its design and its relation 
to interpolation. It is well-known that high-order approximations to a discontinuous 
function generate spurious oscillations near the discontinuity (Gibbs phenomenon). 
Unlike standard finite-difference methods which use a fixed stencil, modern shock- 
capturing schemes use an adaptive stencil which is selected according to the local 
smoothness of the solution. Near discontinuities this technique automatically switches 
to one-sided approximations, thus avoiding the use of discontinuous data which brings 
about spurious oscillations. 
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1 Introduction 


In this paper, we describe and analyze numerical techniques that are designed to 
approximate weak solutions of hyperbolic systems of conservation laws in several 
space dimensions. For sake of exposition, we shall describe these methods as they 
apply to the pure initial value problem (IVP) for a one- dimensional scalar conservation 
law 

«* + /(u)» = 0, u(x, 0) = u 0 (x). (1.1) 

To further simplify our presentation, we assume that the flux f(u) is a convex function, 
i.e., f"(u ) > 0 and that the initial data uo(x) are piecewise smooth functions which 
are either periodic or of compact support. Under these assumptions, no matter how 
smooth uq is, the solution u(x,t) of the IVP (1.1) becomes discontinuous at some 
finite time t = t c . In order to extend the solution for t > f c , we introduce the notion 
of weak solutions, which satisfy 

dt L u dx + *)) ~ /( u ( a > *)) = 0 (1.2a) 

for all 6 > a and t > 0. Relation (1.2a) implies that u(x,t ) satisfies the PDE in (1.1) 
wherever it is smooth, and the Rankine-Hugoniot jump relation 

f(u(y + 0, t)) - f(u(y - 0, t )) = [u(y + 0, t) - u(y - 0, t)] J (1.26) 

across curves x = y(t) of discontinuity. 

It is well-known that weak solutions are not uniquely determined by their initial 
data. To overcome this difficulty, we consider the IVP (1.1) to be the vanishing 
viscosity limit t \ 0 of the parabolic problem 

( u ‘)t + f(u e )x = eu c xx u'(x,0) = uo(x), (1.3a) 

and identify the unique “physically relevant” weak solution of (1.1) by 

u = ]W. (1.36) 

The limit solution (1.3) can be characterized by an inequality that the values ul = 
u (y ~ 0>0i U R ~ u {v + 0) i) a ^d 3 = dy/dt have to satisfy; this inequality is called an 
entropy condition; admissible discontinuities are called shocks. When f(u) is convex, 
this inequality is equivalent to Lax’s shock condition 

o(«l) > s> a(u R ) (1.4) 

where a(u) = f'(u) is the characteristic speed (see [8] for more details). 

We turn now to describe finite difference approximations for the numerical solution 
of the IVP (1.1). Let v" denote the numerical approximation to u(xj, t n ) where 
x j = jh , t n = nr; let Vh(x, t ) be a globally defined numerical approximation associated 
with the discrete values {u"},oo < j < oo,n > 0. 
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The classical approach to the design of numerical methods for partial differential 
equations is to obtain a solvable set of equations for { v ?} by replacing derivatives in 
the PDE by appropriate discrete approximations. Therefore, there is a conceptual 
difficulty in applying classical methods to compute solutions which may become dis- 
continuous. Lax and Wendroff [9] overcame this difficulty by considering numerical 
approximations to the weak formulation (1.2a) rather than to the PDE (1*1)* For this 
purpose, they have introduced the notion of schemes in conservation form: 

*>; +1 = V? - A(7 i+ J - 7,-1 ) = (Eh ■ (l*a) 

here A = rjh and f i+ 1 denotes 

7i+j = /(»?.*« <+*); (i-») 

. . . ,w 2 k) is a numerical flux function which is consistent with the flux f{u), in 

the sense that __ 

f(u,u,...u) = /(u); (1.5c) 

Eh denotes the numerical solution operator. Lax and Wendroff proved that if the 
numerical approximation converges boundedly almost everywhere to some function 
it, then u is a weak solution of (1.1), i.e. , it satisfies the weak formulation (1.2a). 
Consequently discontinuities in the limit solution automatically satisfy the Rankine- 
Hugoniot relation (1.2b). We refer to this methodology as shock-capturing (a phrase 
coined by H. Lomax). 

In the following, we list the numerical flux function of various 3-point schemes 
{k = 1 in (1.5b)): 

(i) The Lax-Friedrichs scheme [7] 


f(wi,w 2 ) = hf(w i) + f(w 2 ) - ~(w 2 -u;i)] 

(1.6) 

(ii) Godunov’s scheme [1] 


J(w u w 2 ) = f{V(0\ w 1 ,w 2 ))] 

(1.7a) 

here V(x/t] u>i, w 2 ) denotes the self-similar solution of the IVP (1.1) 

data , 

. . J wi x < 0 
U °W = { x>0 ' 

with the initial 
(1.76) 

(iii) The Cole-Murman scheme [12]: 


J(wi,w 2 ) = hf(w i) + f(w 2 ) - |o(uji,u;2)|(u 72 -“U/l)] 

(1.8a) 

where / if W X ± w 2 

= at--.' 

(1.86) 
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(iv) The Lax-Wendroff scheme [9]: 

7K> ™i) = ^{/(^ 1 ) + f(w 2 ) - Ao [/(x^a) - /(xwi)]}- (1.9) 

Let i£(t) denote the evolution operator of the exact solution of (1.1) and let Eh 
denote the numerical solution operator defined by the RHS of (1.5a). We say that the 
numerical scheme is r-th order accurate (in a pointwise sense) if its local truncation 
error satisfies 

E(t) ■ u - E h ■ u = 0(h r+1 ) (1.10) 

for all sufficiently smooth u\ here t = 0(/i). If r > 0, we say that the scheme is 
consistent. 

The schemes of Lax-Friedrichs (1.6), Godunov (1.7), and Cole-Murman (1.8) are 
first order accurate; the scheme of Lax-Wendroff (1.9) is second order accurate. 

We remark that the Lax-Wendroff theorem states that if the scheme is convergent, 
then the limit solution satisfies the weak formulation (1.2b); however, it need not be 
the entropy solution of the problem (see [4]). It is easy to see that the schemes of 
Cole-Murman (1.8) and Lax-Wendroff (1.9) admit a stationary “expansion shock” 
(i.e., /(ul) = /( ur ) with a(ui ) < a(u«)) as a steady solution. This problem can be 
easily rectified by adding sufficient numerical dissipation to the scheme (see [11] and 

[3])- 


2 Interpolatory Schemes and Linear Discontinu- 
ities 

Let us consider the constant coefficient case f(u ) = au,a = const, in (1.1), i.e., 

u t + au x = 0, u(x, 0) = u 0 (x), (2.1 a) 

the solution to which is 

u(x, t) = u 0 (x — at). (2-1&) 

In this case the schemes (1.6) - (1.9) take the form 

»? +, = E ClO-Kw = (E h ■ v"),. (2.2) 

l=-K- 

where v = Xa is the CFL number. The coefficients Ct(v) are independent of the 
numerical solution v n \ this makes Eh a linear operator. 

We say that the numerical scheme Eh is (linearly) stable if 

\\(EhT\\ <C for 0 < nr < T, r = 0(h). (2.3a) 

In the constant coefficient case the scheme is stable if and only if it satisfies von 
Neumann’s condition 

K+ 

| ^2 C t (v)i^\ <1 for all 0 < ( < ir. (2.3 b) 

l=-K- 
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It is easy to see that all the 3 point schemes (1.6) - (1.9) are stable under the CFL 
condition 

|i/| = | Aa| < 1. (2.3c) 

The notion of stability (2.3a) is related to convergence through Lax’s equivalence 
theorem, which states that a consistent linear scheme is convergent if and only if it 
is stable (see [13] for more details). 

Let us denote by Sf the stencil of (r + 1) successive points starting with x,- 

*57 *®«+ 1 > ■ ■ ■ > ®»+r} > (2.4a) 


let P(x] S;,u) denote the unique polynomial of degree r interpolating the (r+1) values 
of u on this stencil and let Q{x\ u ) denote the piecewise polynomial interpolation of 
u 

Q(x]u) = P{x ; S^y, u ) Xj_i < x < Xj (2.46) 

We refer to the numerical scheme 


v ; +1 - 


= Q( x j - CLT]V n ) 


(2.4c) 


as interpolatory scheme. Clearly, the interpolatory scheme (2.4) is r-th order accurate. 
When Q(x\v n ) is the piecewise linear interpolation of v n (i.e., r = 1 , i(j ) = j — 1 
in (2.4b)) then (2.4c) is the first-order accurate upwind scheme; in the constant 
coefficient case this scheme is identical to those of Godunov (1.7) and Cole-Murman 
( 1 - 8 ). 

Next let us assume a > 0 and consider the second order case r = 2 in which 
Q(x] v n ) is a piecewise-parabolic interpolation of v n . There are two different choices of 
stencil in (2.4): Taking Q in to be the parabola through S 2 -_ x = {x,_i,x,-, x J+1 } 

(i.e., i(j ) = j - 1) results in the Lax-Wendroff scheme (1.9); taking Q in [xj-i,Xj\ 
to be the parabola through S*_ 2 = {xj- 2 > x j-i, x j} (i.e., (i(j) — j — 2) results in the 
second-order upwind scheme. 

We turn now to consider the application of these schemes to the step function 


H(x) 


H(x) 


0 x < 0 „ _ f 0 j <0 

1 x > 0 ’ Hj ~ { 1 j > 1 ' 


(2.5a) 


For the first order upwind scheme we get that 


f 0 x < 0 

Q{x]H)= i x/h 0 <x <h ; (2.56) 

[l h < x 


for the Lax-Wendroff scheme 




0 

1 -K 1 - k)(2-S) 
1 


x < —h 
—h < x < 0 
0 < x < h ’ 
h < x 


(2.5c) 
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for the second order upwind scheme we get that 




( 0 


1 + 1(1 - 1)(2 - f ) 
l i 


X < 0 

0 < x < h 
h < x <2h 
2 h < x 


(2.5 d) 


We observe that Q in (2.5b) is a monotone function of x; consequently the numerical 
solution by Godunov’s scheme to these data is also monotone. On the other hand 
Q for the second order schemes (2.5c) - (2.5d) is not a monotone function. For the 
Lax-Wendroff scheme Q is negative in — h < x < 0 and has a minimum of -0.125; 
similarly for the second order upwind scheme Q is larger than 1 in h < x < 2h 
with a maximum of 1.125. This observation explains the Gibbs-like phenomenon of 
generating spurious oscillations in calculating discontinuous data with these second 
order schemes. 

We say that the scheme E \ is monotonicity preserving if 


v monotone =*► Eh ■ v monotone. 


( 2 . 6 ) 


Clearly the numerical solution of a monotonicity preserving scheme to initial data of a 
step-function is always monotone and therefore the discontinuity propagates without 
generating spurious oscillations. 

Godunov has shown that the linear scheme (2.2) is monotonicity preserving if and 
only if 

C t {v) > 0, -K_ <1<K+ ; (2.7) 

this implies that a monotonicity-preserving scheme which is linear is necessarily only 
first-order accurate. It took some time to realize the Godunov’s monotonicity the- 
orem does not mean that there are no high-order accurate monotonicity preserving 
schemes; it only means that there are no such linear ones. Hence high-order accurate 
monotonicity-preserving schemes are nonlinear in an essential way. 

The second-order accurate schemes mentioned above are linear because the choice 
of the stencil (2.4) is fixed. Let us consider now a piecewise-quadratic interpolation 
which is made nonlinear by an adaptive selection of the stencil in (2.4b). For the 
interval [xj_i,Xj] let us consider the two stencils Sj_ 2 = {xy_ 2 , xy-i, Xj} and Sj_i = 
{xj_i, Xj, Xj + i}, and select the one in which the interpolant is smoother. If we measure 
the smoothness of u by the second derivative of the corresponding parabola we select 



3-2 
3 ~ 1 


if |^rP(*;^_ a ,u)| < |£ r P(x;5?_ 1) u)| 
otherwise 


(2.8a) 


When we apply this selection of stencil to the step-function H(x) (2.5a) we get that 
for [x_i, x 0 ] we choose the stencil Si 2 = {^- 2 ) *-ii *o} f° r which P(x; 5’i 2 , H ) = 0; for 
the interval [xj, x 2 ] we choose the stencil = {x 1} x 2 , x 3 } for which P(x; S*, H) = 1. 
As is evident from comparing (2.5c) and (2.5d) it does not matter which stencil we 
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assign to [xc^Xi] since both parabolae are monotone there; with (2.8a) we select 5^ 
for [x 0 , ®i]- Thus we get in (2.4) 


f 0 x < 0 




+ o<x<h 

1 h < x 


(2.8 b) 


which is a monotone function of x although it is actually a piecewise-quadratic poly- 
nomial. 

The use of an adaptive stencil is the main idea behind the Essentially Non- 
Oscillatory (ENO) schemes to be described later in this paper. It extends to high 
order of accuracy in a straightforward manner: For r-th order accuracy we consider 
for Xj] the r stencils Sj_ r , Sj_ T+1 , ... , We choose i(j) in (2.4b) to be the 

one which minimizes 



3 Total Variation Stability and TVD Schemes 


An immense body of work has been done to find out whether stability of constant 
coefficient scheme with respect to all “frozen coefficients” associated with the problem, 
implies convergence in the variable coefficient case and in the nonlinear case. 

In the variable coefficient case, where the numerical solution operator is linear 
and Lax’s equivalence theorem holds, it comes out that the stability of the variable 
coefficient scheme depends strongly on the dissipativity of the constant coefficient 
one, i.e., on the particular way it damps the high-frequency components in the Fourier 
representation of the numerical solution. 

In the nonlinear case, under assumptions of sufficient smoothness of the PDE, its 
solution and the functional definition of the numerical scheme, Strang proved that 
linear stability of the first variation of the scheme implies its convergence; we refer 
the reader to [13] for more details. 

In the case of discontinuous solutions of nonlinear problems, linearly stable schemes 
are not necessarily convergent; when such a scheme fails to converge, we refer to this 
case as “nonlinear instability.” The occurrence of a nonlinear instability is usually 
associated with insufficient numerical dissipation which triggers exponential growth 
of the high-frequency components of the numerical solution. 

The following theorem states that a stronger sense of stability, namely uniform 
boundedness of the total variation of the numerical solution, does imply convergence 
to a weak solution. 


Theorem 3.1. Let Vh be a numerical solution of a conservative scheme (1.5). 

(i) u 

TV(v k (-,t))<CTV( I*,) (3.1) 

where TV( ) denotes the total variation in x and C is a constant independent of h 
for 0 < t < T, then any refinement sequence h — > 0 with r = 0(h) has a convergent 
subsequence hj — > 0 that converges in L l ° c to a weak solution of (1.1). 


i 


I 


i 
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(ii) If Vh is consistent with an entropy inequality which implies uniqueness of the 
IVP (1.1), then the scheme is convergent (i.e., all subsequences have the same limit, 
which is the unique entropy solution of the IVP (1.1)). 

We say that the scheme Eh is Total Variation Diminishing (TVD) if 

TV (E h ■ v) < TV(v) (3.2) 


where 


TV(w) = ki+i ~™;l- 


(3.3) 


Clearly TVD schemes satisfy (3.1) with (7 = 1 and therefore are TV stable. 
In [2] we have shown that if the scheme can be written in the form 

= iff + Ct i A j + 1 v n - CT iA_ iv" 

where (7* i satisfy for all j 

J-r j 


e£.>0, c^ + c- H < 1 


(3.4a) 


(3.46) 


then the scheme is TVD; here — vfo - v?. Applying this lemma to the 

general scheme 


't 


we get that if A q satisfies 

then the scheme (3.5) is TVD; here 

S i+J : 


A (7 i+ i-7i-i) 

(3.5a) 

/;+i - 9;+£ A i+p n ) 

(3.56) 

< H+* < 1 

(3.6a) 

i « 

1-1 ■ c 

4 4 

(3.66) 


A ' . 1 
J+5 


This shows that the Cole-Murman scheme (1.8) for which q = |a| is TVD subject to 
the CFL restriction A|a J+ i| < 1. 

Using conditions (3.4b) it is possible to construct TVD schemes which are second- 
order accurate in the Lj-sense (see [2] and [14]). However, TVD schemes are at 
most second-order accurate (see [5]). In order to design higher-order accurate shock 
capturing schemes we introduce the notion of Essentially Non-Oscillatory (ENO) 
schemes. 


4 ENO Schemes 

In this section we describe high-order accurate Godunov-type schemes which are a 
generalization of Godunov’s scheme (1.7) and van Leer’s MUSCL scheme [10]. 
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We start with some notations: Let {Ij} be a partition of the real line; let A(I) 
denote the interval-averaging (or “cell-averaging”) operator 

A(I) ■ w = jyj Jw{y)dy\ (4.1) 

let Wj = A(Ij ) • w and denote w = {w-}. We denote the approximate reconstruction 
of iw(x) from its given cell-averages {uJj} by R(x\w). To be precise, R(x\w) is a 


piecewise-polynomial function of degree (r — 1), which satisfies 

(i) R(x]w) = w(x ) + 0(h T ) wherever w is smooth (4.2a) 

(ii) A(Ij) ■ R(\w) = Wj (conservation). (4.26) 

Finally, we define Godunov-type schemes by 

v”+' = Ail,) ■ E( t) ■ *(•;«■) = (E h ■ v"), (4.3a) 

v° = A(I,)u o; ( 4 . 36 ) 

here E(t) is the evolution operator of (1.1). 


In the scalar case, both the cell-averaging operator A(Ij) and the solution operator 
E(t) are order- preserving, and consequently also total-variation diminishing (TVD); 
hence 

TV (Eh ■ w) < TV(R(;w)). (4.4) 

This shows that the total variation of the numerical solution of Godunov-type 
schemes is dominated by that of the reconstruction step. 

We turn now to describe the recently developed essentially non-oscillatory (ENO) 
schemes of [5, 6], which can be made accurate to any finite order r. These are 
Godunov- type schemes (4.3) in which the reconstruction R(x\w), in addition to re- 
lations (4.2), also satisfies 

TV(R(-,W)) < TV(w) + 0(h 1+p ), p > 0 (4.5) 

for any piecewise-smooth function u;(x). Such a reconstruction is essentially non- 
oscillatory in the sense that it may not have a Gibbs-like phenomenon at jump- 
discontinuities of w(x), which involves the generation of 0(1) spurious oscillations 
(that are proportional to the size of the jump); it can, however, have small spurious 
oscillations which are produced in the smooth part of w[x), and are usually of the 
size 0{h T ) of the reconstruction error (4.2a). 

When we use an essentially non-oscillatory reconstruction in a Godunov-type 
scheme, it follows form (4.4) and (4.5) that the resulting scheme (4.3) is likewise 
essentially non-oscillatory (ENO) in the sense that for all piecewise-smooth function 
w(x) 

TV (E h ■ w) < TV(w) + 0(h 1+p ), p> 0; (4.6) 

i.e., it is “almost TVD.” Property (4.6) makes it reasonable to believe that the total 
variation of the numerical solution is uniformly bounded. We recall that by Theorem 
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3.1, this would imply that the scheme is convergent (at least in the sense of having 
convergent subsequences). This hope is supported by a very large number of numerical 
experiments. 

Next we describe one of the techniques to obtain an ENO reconstruction. To 
simplify our presentation we assume that {iy} is a uniform partition 

Ij = Xj = jh. 

Given cell averages {toy} of piecewise- smooth function w(x), we observe that 

hWj = f 3 w(y)dy = W(xj) - W(xy_ x ) (4.7a) 

Jxj-l 

where 

W(x ) = / w{y)dy (4.7 b) 

Jxo 

is the primitive function of w(x). Hence we can easily compute the point values 
{W(xy)} by summation 

W(xi) = hJ2wj. (4.7 c) 

3=i o 

Once we have computed the point values of the primitive function we use the ENO in- 
terpolation technique (2.4), (2.9) to obtain Q(x ; W ), an r-th order piecewise- 
polynomial interpolation of W, i.e., 


Q(x] W) = P(x; Sy ( yj, W) for xy_! < x < xy (4.8a) 

where P(x; S* , W) is the unique r-th degree polynomial which interpolates W over 
the stencil S\ = {x,, x, + j, . . . , Xi +r }, and i(j ) is chosen so that 


1^' 


(4.8 b) 

We define R(x\ w) by 

R(x,W) = i-QO; W). 

(4.9) 


We observe that if w(x) is smooth in (xy_ x ,xy) then for h sufficiently small the 
algorithm (4.8b) will select a stencil SJj-y» in which w(x) is smooth. It follows then 
from standard interpolation theorems that 


i2(x; = ^ P(x; %)> W )\ = ^ w + °( hT ) = w ( x ) + °( hr ) ( 4 -! 0 ) 

which is property (4.2a). Furthermore (4.10) holds in every interval except for those 
in which w(x) has a discontinuity. As we have seen in the examples (2.5) and (2.8b) 
the Gibbs-phenomenon is associated with intervals near the discontinuity and not 
with the interval that contains the discontinuity. This is why the reconstruction (4.8) 
- (4.9) satisfies the ENO property (4.5); in [2] we show that the second-order accurate 
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ENO scheme is actually TVD. The conservation property (4.2b) follows directly from 
the definition (4.9): 


A(Ij)R(-‘, w ) 


i [“ -ftyx | W)dx = W) - W)\ 

it */xj_ i ctx ft 

jpyfe) -«"(*<->)] =w<- 


(4.11) 


The abstract scheme (4.3) can be written in the standard conservation form (1.5). 
To do so let us denote by v(x,t) the solution in the small of the IVP 


f (v)t + f(v)x = 0 
( v(x, 0) = R(x; v n ) 


0<<<T 


(4.12) 


and integrate this PDE over Ij x [0,r]; using the divergence theorem and (4.2b) we 
get that u n+1 in (4.3) can be expressed by 

”" + ' = •? - - !j-i\ («•!*•) 

where 

f)+$ = ~j 0 /( s ( x i»0)* (4.136) 

In the first-order case the scheme (4.13) is identical to Godunov’s scheme and the 
numerical flux (4.13b) can be expressed in a closed form by (1.7b). For higher order 
schemes we use a numerical flux which is an appropriate approximation to (4.13b) 
(see [6] for more details). 

We remark that the ENO schemes are related to the interpolatory schemes of Sect. 
2 as follows: In the constant coefficient case a fixed choice of stencil (i.e. , i(j) — j 
= constant in (4.8a)) results in the interpolatory scheme (2.4) corresponding to the 
same choice of stencil. 
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